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Abstract 

We investigate local update algorithms for the fully frustrated XY 
model on a square lattice. In addition to the standard updating pro- 
cedures like the Metropolis or heat bath algorithm we include overre- 
laxation sweeps, implemented through single spin updates that pre- 
serve the energy of the configuration. The dynamical critical exponent 
(of order two) stays more or less unchanged. However, the integrated 
autocorrelation times of the algorithm can be significantly reduced. 
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Figure 1: A ground state configuration of the FFXY model. The 
dotted lines represent links with antiferromagnetic couplings, the full 
lines are links with ferromagnetic couplings. 



1 Introduction 

The 2-dimensional fully frustrated XY (FFXY) model has attracted a lot of 
attention over the years. It is defined through the partition function 

(1) 

The 2-component spins Sj are defined on a square grid. They are constrained 
to unit length, = 1. (3 > denotes the inverse temperature. Units are 
chosen such that ks = 1- The are either +1 or —1 in such a way that each 
elementary plaquette of the lattice contains exactly one antiferromagnetic 
bond (with J = —1). We stick to the convention that the antiferromagnetic 
links are those lying in every second row of the lattice, see figure [I], which 
shows also one of the ground states of the FFXY model. 

The importance of the frustrated XY model stems from the fact that it 
provides a convenient framework to study a variety of interesting phenom- 
ena displayed by numerous physical systems. Experimental realizations of 
the model are, e.g., 2-dimensional arrays of Josephson junctions or supercon- 
ducting wire networks JTJ, ^ [|. 

The phase structure of the FFXY model is still debated. In the literature, 
there seems to be agreement about the following: In addition to the type of 
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excitations present in the unfrustrated XY model [|j] (spin waves at large ft 
and vortices for ft < (3kt) there is a new type of discrete excitations which 
should be relevant at sufficiently large ft. With each elementary plaquette i 
one can associate a chirality variable Oi, defined through 

Oi = sign (sin (f) 12 + sin 2 3 + sin 34 + sin 4 i) , (2) 
0fcz = 9 k — 61 — A H . 

The corners of the plaquette i are labelled clockwise 1 ... 4. The angles 
8k are defined through the relation Sk = (cos 9k, siia 9k), and Am = or n 
for ferromagnetic and antiferromagnetic links, respectively. <Tj can assume 
values ±1. The global symmetry of the model with respect to the reflection 
of all a's can be spontaneously broken. If this transition takes place at a ft- 
value different from that where the XY degrees of freedom undergo a (most 
likely) Kosterlitz-Thouless phase transition, then it follows from standard 
wisdom about universality that it should belong to the universality class of 
the 2-dimensional Ising model. On the other hand, if the two transitions 
are exactly on top of each other, one expects a new universality class, with 
exponents different from the 2D Ising values. 

The critical properties of the FFXY model have been investigated in a 
number of Monte Carlo studies, see, e.g., refs. p||-fl6|. All studies, including 
the most recent ones, indicate that there is an urgent need to go to larger 
lattices. Unfortunately, the simulations of the FFXY model exhibit strong 
critical slowing down. On large lattices one has to perform many update 
sweeps in order to obtain a statistically independent configuration. E.g., at 
ft = 2.2, which is in the critical region, the autocorrelation time of a local 
Metropolis algorithm is ~ 1600 on a 64 by 64 lattice. Fortunately, there seems 
to be no exponential slowing down. At criticality the autocorrelations grow 
like t ~ AL Z , where L denotes the linear lattice extension. The dynamical 
critical exponent is about 2.3 for a local Metropolis update algorithm. 

It is tempting to apply cluster algorithms [17|, [18|] to the problem. How- 
ever, as has been experienced quite often in the case of frustrated models, 
there is no obvious way to avoid cluster size distributions with strong weight 
on very large cluster sizes, thus spoiling any gain in efficiency. Using the 
general formalism of Kandel et al. ||19|| , with the useful extensions of ref. |20 



one can convince oneself that (distinct from the case of the fully frustrated 
Ising model |21]]) in a quite large class of algorithms there is none that avoids 
this problem. 
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In ref. || an attempt was made to improve on the performance of the 
usual local heat bath or Metropolis algorithm, employing Fourier accelera- 
tion. However, there was negligible speedup of observables built from the 
discrete a-variables introduced above. In ref. a Hybrid Monte Carlo algo- 
rithm was used, but without a quantification of critical slowing down. 

We here investigate the inclusion of microcanonical overrelaxation steps 
p2|-p9| in the update procedure. To the best of our knowledge it has not yet 
been employed in the context of the frustrated XY model. The algorithm is 
very easy to implement, keeps the updating perfectly local and thus allows 
simple vectorization and/or parallelization. Still, an effective speedup of one 
order of magnitude can be obtained. 



2 Observables and Algorithms 

In order to test a number of algorithms, we measured for both the XY degrees 
of freedom s and the Ising type variables a the energies, 

Ei = (jMfij)- (3) 

Here, i and j are nearest neighbours in the lattice, and the variables fii are 
defined through 

in = (-iy^y Gi , (4) 

where the lattice coordinate i has components (i x , i y ). Furthermore, we mea- 
sured the susceptibilities 

Xxy = ^2(8i • Sj) , 
3 

Xi = ^2(fhlij)- ( 5 ) 

j 

% is any lattice site, and j runs over the whole grid. All our simulations were 
done on L by L lattices, L even, with periodic boundary conditions. 

In order to quantify the critical slowing down we measured the integrated 
autocorrelation times, defined through 

J oo 

Tint,X = ~Y. C Xx(t) , (6) 
1 t=0 
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obs 


Tint 


Et 


9.96(16) 


Xi 


11.16(20) 




8.36(12) 


XXY 


8.12(12) 



Table 1: Integrated autocorrelation times on an 8 x 8 lattice, heat 
bath algorithm, (3 = 2.2. 



using the selfconsistent window method proposed in |3(|. Cxx denotes the 
normalized autocorrelation function of the observable X. 

Table [l] shows the integrated correlation times of the four observables 
defined above, for L = 8. Here and in the following we chose (3 = 2.2. This 
value is compatible with most estimates of the critical inverse temperature in 
the literature. Our results for r int were obtained from 10 6 lattice sweeps with 
a heat bath algorithm, xi has the longest autocorrelations. This remains 
so also on other lattice sizes and also with any other algorithm that we 
investigated. 

In the following, we shall report on results obtained with four types of 
algorithms. In all cases the sweeps are done in a typewriter fashion. 

Metropolis (M): 

At site i, a rotation of the spin by an angle 9 is proposed, where 9 is cho- 
sen with uniform distribution from the interval [— e, +e]. e can be used to 
tune the acceptance rate. A Metropolis hit is accepted with probability 
min[l, exp(— /3AE)], where AE denotes the change of energy associated with 
the proposed rotation. At site i, one performs n^t Metropolis hits. 

Metropolis Reflection (MR): 

This is a Metropolis algorithm based on the operation which lies at the 
heart of the Wolff cluster algorithm for the XY model. A single random 
unit vector u is selected. One sweeps through the lattice. At site i, the 
spin Si is decomposed in components parallel and perpendicular to u. It 
is proposed to change the sign of the parallel component. The reflection is 
actually performed with a probability determined by the Metropolis rule. 

Heat Bath (HB): 

We employ a heat bath routine for the generation of U(l) random numbers 
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developed by Hattori and Nakajima [51], based on a rejection method with 
transformation of variables. It has a high acceptance rate for all values of 
the temperature. 



Micro canonical Overrelaxation (OR): 

Single spins are moved such that the energy remains unchanged. The con- 
tribution to the total energy of the system that depends on a spin si is 
proportional to s~l ■ hi, where hi is a weighted sum of the nearest neighbour 
spins of si. The spin s~l is decomposed into components parallel and per- 
pendicular to hi. The latter component is reflected (multiplied by —1) with 
probability 1, leaving unchanged the scalar product s~l • hi. Note that the OR 
algorithm is non-ergodic. It can only be used in a mixture with an ergodic 
algorithm. 

3 Metropolis and Heat Bath Algorithms 

We started our investigation with a study of the first three algorithms to 
compare their efficiency. In table ^] we quote autocorrelation times for the 
"slowest" observable Xi f° r the Metropolis- Algorithm with various choices of 
parameters, for the Metropolis reflection and for the heat bath algorithm. 
The simulations were done at (3 — 2.2, on an 8 by 8 lattice. We always 
performed one million full lattice sweeps. In the Metropolis cases we also 
quote acceptance rates. In all cases we measured the CPU time consumed 
for 20,000 sweeps on a Pentium 166 MMX PC. The last column contains 
a number r e g- which is obtained by multiplying the autocorrelation time by 
the CPU factor, and then normalizing such that r e g- is one for the heat bath 
algorithm. 

The heat bath algorithm has by far the smallest autocorrelation times. 
However, the Metropolis reflection algorithm implementation is so much 
faster that it effectively becomes the most efficient algorithm among the ones 
investigated. In the combination with the OR algorithm to be discussed in 
the next section we always used two subsequent Metropolis reflection algo- 
rithm sweeps as the basic unit. Let us remark that using different algorithm 
implementations or computers we could have reached different conclusions. 

We now investigated the critical slowing down of the Metropolis reflection 
algorithm by running it on various lattice sizes, always at (5 = 2.2. Table ^ 
shows our findings for lattice sizes up to L = 64. The integrated autocorrela- 
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algorithm 


e 


nut 


T int,Xl 


acc 


cpu 




M 


1.5 


1 


43.00 ± 1.40 


47% 


7.1 


0.943(31) 


M 


2.5 


1 


34.40 ± 1.00 


30% 


7.3 


0.776(23) 


M 


3.0 


1 


35.92 ± 1.04 


26% 


7.4 


0.821(24) 


M 


3.0 


2 


20.40 ± 0.44 


51% 


14.2 


0.895(19) 


M 


3.14 


1 


36.00 ± 1.08 


25% 


7.4 


0.823(25) 


M 


3.14 


2 


21.56 ± 0.48 


48% 


14.5 


0.966(22) 


MR 






29.68 ± 0.56 


25% 


4.1 


0.376(07) 


HB 






11.16 ± 0.20 




29.0 


1.000 



Table 2: Integrated autocorrelation times of Xi> L = 8. acc denotes 
the acceptance rate, cpu is a rough estimate of CPU, in seconds, 
for 20,000 sweeps. r e g- is obtained from multiplying the integrated 
autocorrelation time by the cpu factor and normalizing such that it is 
one for the heat bath algorithm. 



L 


Tint,_B XY / L 2 ' 3 


7int,xi/-k 2 ' 3 


S/10 y 


8 


0.1051(17) 


0.1306(21) 


1.3 


12 


0.0763(17) 


0.1082(28) 


1.3 


16 


0.0760(17) 


0.1092(27) 


2.5 


20 


0.0753(17) 


0.1050(27) 


3.8 


24 


0.0813(17) 


0.1112(28) 


6.3 


32 


0.0825(17) 


0.1074(31) 


9.0 


40 


0.0798(17) 


0.1080(38) 


10.0 


48 


0.0796(17) 


0.1087(43) 


12.0 


64 


0.0871(17) 


0.1117(60) 


13.7 



Table 3: Integrated autocorrelation times of E XY and \i f° r the 
Metropolis reflection algorithm, rescaled by a factor L 2 3 in both cases. 
The correlation times are measured in units of two sweeps. S denotes 
the statistics in the same units. 
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tion times have been rescaled by a factor L 2 3 . The T int become really large. 
For L = 64 we observe t m t y, = 1593 ± 86. Fits with the law 

r int = A ■ U (7) 

yield 

A z x 2 / d °f 

E XY 0.064(4) 2.36(2) O 
Xi 0.106(8) 2.31(2) 0.5 

2 Metropolis Reflections 

In both cases the data from L = 8 were discarded from the fit. Our error 
bars take into account statistical errors only. 

4 Efficiency of Overrelaxation Steps 

We decided to examine a combination of two Metropolis reflection sweeps 
with iV OR sweeps. It is a general observation that N should scale like the 
correlation length or, if the system is at criticality, with the lattice size L. 
The latter case is relevant for us. 

In order to find the right mixture of reflection and OR steps, we computed 
ri n t for the XY energy and the Ising susceptibility. The runs were done on 
SUN spare stations, for lattice sizes 8, 16, and 32. Our results for the two 
smaller lattices are quoted in tables |4] and [|. In the particular implementation 
of our code, the overrelaxation sweeps needed m 1/16 of the run time of two 
MR sweeps. This CPU factor enters into the numbers t^°" defined through 

C" = 71* (l + ^) • (8) 

The behavior of rf°™ in both tables together with the L = 32 results suggest 
that it is reasonable to mix according to N — | L . Of course this would have 
come out differently for other CPU time ratios between the MR and the OR 
parts of the update procedure. On a Pentium PC, we measured a factor of 
1/6 between the two CPU times. 

In order to study the volume dependence of the r int we made > 500, 000 
update cycles on various lattice sizes, ranging from L = 8 to L = 128. A 



7 



N 


7int,i?xY 


_corr 
'int,_B X Y 


T int,Xl 


,_corr 
int,xi 





9.54±0.26 


9.54±0.26 


11.88±0.41 


11.88±0.41 


1 


5.82±0.13 


6.18±0.14 


5.63±0.12 


5.98±0.13 


2 


5.50±0.12 


6.19±0.14 


4.84±0.09 


5.45±0.10 


3 


5.39±0.11 


6.40±0.13 


4.50±0.08 


5.34±0.10 


4 


5.39±0.11 


6.74±0.14 


4.41±0.08 


5.51±0.10 


5 


5.18±0.11 


6.80±0.14 


4.21±0.07 


5.53±0.09 


6 


4.88±0.09 


6.71±0.12 


3.94±0.08 


5.42±0.11 


7 


4.79±0.09 


6.89±0.13 


3.81±0.06 


5.48±0.09 


8 


4.76±0.09 


7.14±0.14 


3.68±0.06 


5.52±0.09 



Table 4: Integrated autocorrelation times for E XY and L = 8. 
Two MR sweeps are blended with N OR sweeps, f3 = 2.2. 



JV 


"HntjExY 


corr 
'int.ExY 


T int,xi 


corr 
int,xi 





32.70±1.87 


32.70±1.87 


56.94±3.95 


56.94±3.95 


1 


13.25±0.43 


14.08±0.46 


20.29±0.82 


21.56±0.87 


2 


10.57±0.31 


11.89±0.35 


13.89±0.46 


15.63±0.52 


3 


9.60±0.26 


11.40±0.31 


10.85±0.32 


12.88±0.38 


4 


9.08±0.26 


11.35±0.33 


10.24±0.30 


12.80±0.38 


5 


8.60±0.23 


11.29±0.30 


8.96±0.24 


11.76±0.32 


6 


7.79±0.19 


10.71±0.26 


7.87±0.20 


10.82±0.28 


7 


7.72±0.23 


11.19±0.33 


7.61±0.21 


10.94±0.30 


8 


7.67±0.21 


11.51±0.32 


7.38±0.17 


11.07±0.26 


9 


7.69±0.19 


12.02±0.30 


7.20±0.17 


11.25±0.27 




Table 5: 


Same as table 


H but for L ■- 


= 16. 
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L 


Tmt,E XY / L 




S/10 y 


8 


0.514(10) 


0.2717(48) 


0.5 


12 


0.430(10) 


0.2343(56) 


0.5 


16 


0.340(08) 


0.1864(47) 


0.5 


20 


0.341(10) 


0.1926(58) 


0.5 


24 


0.322(10) 


0.1822(60) 


0.5 


32 


0.289(12) 


0.1656(62) 


0.5 


40 


0.296(12) 


0.1694(76) 


0.5 


48 


0.290(13) 


0.1714(86) 


0.5 


64 


0.262(14) 


0.1564(98) 


0.5 


80 


0.291(08) 


0.1720(57) 


2.5 


96 


0.297(09) 


0.1727(62) 


2.5 


128 


0.288(10) 


0.1606(69) 


1.5 



Table 6: Integrated autocorrelation times for E^ Y and Xi °f the 
Metropolis reflection algorithm, blended with N = 3L/8 OR sweeps. 
The autocorrelation times have been rescaled by factors L 113 and 
L 1-35 , respectively. A unit of measurement is two Metropolis reflection 
sweeps combined with 3L/8 OR sweeps. S denotes the statistics in 
these units. 
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cycle consisted of two MR sweeps, followed by 3L/8 OR sweeps. Our results 
for the r int are quoted in table |[ 

Fitting with the law eq. (|7|), discarding lattice sizes L < 24, we obtained 

A z X 2 / d °f 
E XY 0.28(5) 1.13(4) L0 
Xi 0.17(3) 1.34(4) 0.8 

2 Metropolis Reflections + 3L/8 OR 



For a fair comparison with the Metropolis reflection algorithm one needs 
to use units that take into account the work spent on the OR steps. In fair 
work units the scaling law is 

r int = A-U (l + u;3L/8) , (9) 

with 

1 CPU cost of OR sweep . . 

U= 2 CPU cost of MR sweep ' ^ ' 

The asymptotic behaviour will thus be governed by an exponent 2+1, with an 
amplitude 3/8 to A. Comparing the results of the OR and the pure Metropo- 
lis algorithms, we find compatible or similar exponents. The effective ampli- 
tudes, however, differ. In case of the chiral susceptibility, the net efficiency 
gain for large lattices is 

3 0.17(3) 

7 = - U) 7+. 11 

' 8 0.106(8) V ; 

For realistic values of uj between 0.1 and 0.2 this gives a speedup of about 
10. This statement is in good agreement with the actual values of autocor- 
relations times observed in our simulations. 

It is interesting to look also at the performance of an algorithm with a 
fixed admixture of OR steps. We made simulations with N = 3 kept fixed 
for all lattice sizes. Our results are summarized in table [7|. 

Fitting with the scaling law eq. (0), we had to discard data from L < 32. 
The results are 

a z rV d °f 

Exy 0.10(1) 2.33(3) 22 
Xx 0.022(3) 2.21(3) 2.5 

2 Metropolis Reflections + 3 OR 
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L 


Tint Evv/i 2 ' 32 


* lilt Yl / 


S/10 e 


8 


0.04169(96) 


0.04391(93) 


0.3 


16 


0.01679(39) 


0.02867(72) 


0.5 


20 


0.01495(24) 


0.02788(55) 


1.3 


24 


0.01419(19) 


0.02722(46) 


2.5 


32 


0.01066(19) 


0.02280(49) 


2.5 


36 


0.00988(20) 


0.02235(53) 


2.5 


40 


0.01004(23) 


0.02388(66) 


2.5 


48 


0.01029(20) 


0.02385(56) 


5.0 


64 


0.00945(35) 


0.02080(90) 


2.5 


72 


0.01048(27) 


0.02362(70) 


8.0 


96 


0.00961(42) 


0.02315(18) 


5.0 



Table 7: Integrated autocorrelation times for E XY and Xi f° r the 
Metropolis reflection algorithm, blended with N = 3 OR sweeps. The 
autocorrelation times have been rescaled by factors L 2 23 and L 2 20 , 
respectively. A unit of measurement is two Metropolis reflection steps 
combined with 3 OR steps. S denotes the statistics in these units. 



Assuming an u about 0.15, the net efficiency gain over the Metropolis 
reflection without OR steps is of order 3 (for Xi)- 

Conclusion 

The net efficiency gain of an inclusion of OR steps in a Monte Carlo proce- 
dure for the FFXY model might vary with implementation and fine tuning 
details. It does not seem to improve on the dynamical critical exponents 
in a significant way. However, for a given lattice size, it definitely enhances 
the efficiency of the updating. The fact that it can be implemented so eas- 
ily makes its usage mandatory. It would be interesting to combine the OR 
updating with multilevel techniques. E.g., one could perform simultaneous 
rotations of all spins in a block of size I by I such that the total energy 
remains unchanged. 
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